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SUMMARY 


Sand production is a complex physical process that depends on the external stress and flow rate conditions as 
well as on the state of the material. Models developed for the prediction of sand production are usually 
solved numerically because of the complexity of the governing equations. Testing of new sand production 
models can very well be performed through calibration with laboratory experiments, which by construction 
possess geometric symmetry facilitating explicit mathematical analysis. We introduce an erosion model that 
is built upon the physics (poro-mechanical coupling of the fluid-solid system) usually incorporated in ero- 
sion models for the prediction of sand production. Around this model, we set up a mathematical framework 
in which sand production models because of erosion can be tested and calibrated without having to resort to 
complex numerical work or specialised software. The model is validated by data of volumetric sand produc- 
tion from a hollow cylinder test on synthetic sandstone. Generalisations of the model, which are naturally 
incorporated in the same framework and have useful phenomenological features, are discussed. Copyright 
© 2015 John Wiley & Sons, Ltd. 
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1. INTRODUCTION 


Sand erosion, in the sense that interests us in this work, is the phenomenon where a fluid saturated rock 
loses its mechanical integrity because of the applied stress field and is transported away in the presence 
of pore pressure gradients. The unintended by-product of solid particles of oil and gas production is 
generally referred to as sand production. In practice, the aggressive pumping of oil from the well 
causes grain dislocation from the solid matrix of the rock, thereby leading to mechanical problems 
such as accumulation of sand in the wellbore and the creation of unstable cavities in the geological 
formation. Particle influx into the wellbore may lead to various problems such as erosion of valves 
and pipelines, plugging of production liners and sand deposits in the separators. Additionally, during 
production, sudden erosion in high-pressure gas wells represents a major safety risk [1, 2]. The 
design issue that arises from this problem, and it is of considerable interest to the oil and gas 
operating companies, is the prediction or control of these particles. 

In the view of modelling, the different processes that are involved in the sand production problem 
are associated with (1) fluid and solids transport; (2) fluid/rock interaction; and (3) rock deformation. 
Modelling of the sand production physical problem is not a trivial process as it involves the 
following coupled mechanisms: (a) the mechanical instabilities and localised compressional or 
tensile failure (damage) of the rock in the vicinity of the wellbore because of stress concentration 
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and (b) the hydro-mechanical instabilities caused from internal and surface erosion, which manifest 
themselves in releasing and transferring the particles caused by the action of seepage forces. These 
two mechanisms are coupled because the stress concentration will lead to loss of cohesion of the 
material causing localised damage, which in turn will increase the amount of loose solids that can 
leave the rock formation. This solid wash-out will bring an increase at the porosity of the intact rock 
as a direct effect of the localised damage and eventually will cause the inter-granular forces to re- 
adjust to the new configuration of the rock formation, thus leading to further damage of the rock [3]. 

In reality, erosion is influenced by various mechanisms. These include combinations of drawdown, 
drawdown rate (ramp-up strategy), depletion, flow rate, water-cut, completion strategy (size, phasing 
and orientation of perforations) and frequency of shut-downs and start-ups. It is evident that modelling 
the contribution from all these mechanisms in a single erosion model is highly complex [4—6]. 

It has been suggested that some limited amount of sand production can eventually lead to an 
increase in oil and gas production [7]. This suggestion dictates the necessity to determine the 
amount of sand produced after the first occurrence. This has motivated many researchers to focus 
on the sand prediction via numerical, analytical and experimental investigations. The availability 
of models to predict the onset and quantity of sand produced is very crucial. There is a vast 
literature in the modelling of erosion phenomena attempting to capture the failure and onset of 
sanding [8-19]. Sanding criteria that are routinely used in the models for predictions are usually 
based on shear, tensile failure [9, 10] critical pressure gradient [11], critical plastic deformation 
[12, 13] and erosion based criteria [14-16]. Among those contributions, the authors of 
Papanastasiou and Vardoulakis [20] used bifurcation theory to simulate the localised mechanical 
response of the reservoir rock with evolving pressure and effective stresses during the life span 
of the wellbore into the post bifurcation regime. They showed that the tendency for compressive 
or tensile failure around a cavity depends on the cavity size. Large cavities (e.g., wellbores) fail 
under compression at low stress level, whereas small cavities (e.g., perforations) fail under both, 
extension or compression at higher stress and strain level depending on the material properties. 

The seminal work of Vardoulakis et al. [3] proposed the basic theory for hydrodynamic erosion to deal 
with the instabilities caused from internal and surface seepage forces on sandstone, which is based on 
filtration theory without solving the stress equilibrium equation. In their work, stresses and rock 
deformation were suppressed with emphasis put on mass transport. As an extension to that work, the 
authors of Papamichos and Stavropoulou [12] combined the evolution of localised deformation with 
hydrodynamic erosion calibrated with laboratory experiments to predict the onset and amount of solids 
produced as well as the production rate. Furthermore, the authors of Stavropoulou et al. [21] extended 
the approach by solving the momentum and continuity equations for fluid flow, solids and fluidised 
solids coupled with an elastic perfect plastic Mohr—Coulomb yield criterion to perform wellbore stability 
analysis after erosion. The parameters of cohesion, Young modulus, permeability and sand production 
coefficient were linked with the porosity through a set of calibration parameters. That was the motivation 
for many researchers to adopt the full strength behaviour of the material in their models [13, 22-24]. 

In general, standard analytical models [17, 25] provide formulations for the fluid flow rate that is 
required to cause compressive or tensile failure triggering the onset of sand production. The critical 
fluid flow rate, that onsets the phenomenon, is a function of (a) the viscosity of the permeating fluid, 
(b) the permeability of the rock formation and (c) the formation strength properties. Most numerical 
models [5, 26] deal with the coupling of the physical processes by attempting to link failure to 
plastic yielding through the assumption that the post peak strength is an appropriate condition for 
sand production and tie solids production to a mobilised critical strain level, which is determined by 
laboratory or field data [27, 28]. The use of this criterion requires a calibration study of every 
specific field or laboratory case. However, this criterion based on failed material makes progression 
of failure possible and release more material for production. A detailed general literature survey on 
sand prediction models can be found in Rahmati et al. [29]. 

This paper presents the development of a mathematical framework, in which hydro-mechanical 
models used for predicting sand production in weak rocks (e.g. sandstone) can be tested and 
calibrated. Simultaneously, a simple erosion model is introduced to operate in this framework. The 
model is a transcription of the model given by Papamichos et al. [13]. The hydro-mechanical 
models around our mathematics are constructed along the following line of thought. 
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One starts with a basic phenomenological relation, which states that the rate of change of porosity 
throughout the material depends on the porosity itself and the Darcy flux of the permeating fluid [3]. 
This relation contains a phenomenological coefficient, which we shall call as the erosion strength 
coefficient 4. When 4 is a constant, the basic erosion model equation and the associated models are 
completely hydrodynamic in nature [3, 21]. It turns out that these models cannot localise the effects 
of erosion near the exit surface of the flow in contrast with intuition and observations. A remedy to 
this problem was first suggested in Papamichos and Stavropoulou and Papamichos et al. [12, 13]. It 
is a priori expected that erosion will be stronger in high stress concentration regions. More 
concretely, one may conveniently associate erosion with plasticity. Adopting a certain failure 
criterion for the material under consideration, erosion is restricted in the region where this criterion 
is satisfied; that is, the erosion region is associated with the plastic region [12, 13]. This is 
performed by assuming that the erosion strength depends on the equivalent plastic strain in a 
prescribed manner. The idea is that erosive failure rests on plastic failure. Under plasticity, the 
material is damaged and decohesioned so that is can be eroded away under the weak hydrodynamic 
forces. The procedure of Papamichos and Stavropoulou and Papamichos et al. [12, 13] does indeed 
lead to a working method and reasonable localisation of the erosion effects near the exit surface. 

The same physical ideas can be implemented in a mathematical description of the problem in a 
variety of ways. In this work, the localisation effect is introduced directly in physical space by 
introducing in the model equation a profile factor (a smoothed step function), which depends on the 
size of the plastic region. In other words, the original erosion strength 2 of Vardoulakis et al. [3] 
acquires explicit space dependence so that this parameter ‘knows’ where the plastic region is. This 
is a very convenient simplification of the localisation method proposed in Papamichos and 
Stavropoulou and Papamichos et al. [12, 13], and at least equally as effective. Also, this is how the 
mechanics of the material is introduced in the present erosion model. Degradation of the material, 
which is also present in the previous versions of the erosion model, is formulated suitably to 
complete the coupling of the hydrodynamics with the mechanics of the material. As porosity 
increases in the plastic/erosion region, the material softens (degradation) and the plastic/erosion 
region progressively evolves deeper. The system evolves in a quasi-steady state. The material 
evolves through a succession of states where the equilibrium equations are applied at each step, 
which allows erosion poro-elastoplastic solutions for the stress field to be derived under the assumed 
(cylindrical) symmetry. Thus, a mathematical framework may indeed be set up, in which one can 
test and calibrate various erosion models without having to resort to far more complex numerical 
work (e.g. finite elements) and the use of specialised software. 

The problem is formulated based on the hollow cylinder test geometry. The formulation of the 
differential equations of the erosion kinematics, fluid flow and stress equilibrium is presented in 
Section 2. Comparison of the model predictions with the experimental results of Papamichos et al. 
[13] is presented and discussed in Section 3. The influence of the model parameters on the erosion 
process are discussed in Section 4. In Section 5, the erosion strength 1 is deduced from the 
experimental data by means of back analysis, and the findings are discussed and evaluated. In 
section 6, we comment on certain deficiencies of the basic erosion model and the need for a specific 
type of modifications required. We then introduce suitable generalised models showing that they can 
be naturally included in the set-up mathematical framework. A discussion on the general 
conclusions drawn from this work is given in Section 7. 


2. COUPLED STRESS—FLUID FLOW EROSION MODEL 


The coupled stress—fluid flow erosion model is developed in order to analyse the hollow cylinder test, 
as mentioned earlier. Figure 1 shows the 2D representation of the model. The model is cylindrically 
shaped with a coaxial cylindrical cavity. The fluid flows radially from the outer surface towards to hole. 

Under such conditions, any deformation of the rock will take place in a plane normal to the cylinder 
axis. The model is loaded isotropically at the outer and inner boundaries at rin and Fout With cin and Gout 
respectively. 
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Figure 1. Schematic of the hollow cylinder problem to simulate the erosion process under radial flow. 


2.1. Darcy’s law and fluid mass conservation 


Consider radial flow between two fixed radii, rin <r</rout, in a porous material of thickness H. Fluid 
mass conservation is expressed by the continuity equation: 


ag; 


ax; a 


Darcy’s law for the flow in the porous medium reads 


Op u 

Ox; k( h) qi ( ) 
where q; is defined as minus the Darcy’s velocity for convenience; u is the dynamic viscosity of the 
fluid and k is the intrinsic permeability of the rock, which is a function of porosity @ according to 
the Kozeny—Carman equation: 


p 


k =k, ———_j 
(1-4) 


(3) 


k, is the Kozeny—Carman permeability parameter, which is constant for a given rock. Because of the 
radial symmetry, the mass conservation (Eqn (1)) becomes 


10 
par = 4 
ror (ra) 4) 
where q is radial fluid flux, with the solution 
A(t 
949 ©) 


for some function A(t). The function A(®) is related to the flow rate Q(t) by the relation Q(t) =27HA(t). 
The relation of A(t) with the pressure drop Ap = Pout — Pin is determined by substituting Eqn (5) in the 
Darcy’s law and integrating between the inner and outer radii: 
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=a GE (6) 


2.2. The erosion model 


In a hydro-mechanical context, the rate of erosion is expected to depend on the porosity ¢ of the rock, 
the concentration of the fluidised particles and the Darcy flux q of the fluid flow, or equivalently 
through the Darcy law, the pressure gradient. This dependence can be set in a variety of ways. In 
the original reference [3], the eroded mass rate law is assumed to be proportional to the solid mass 
proportion 1—@, the fluidised particles concentration and the Darcy flux; the proportionality 
constant has dimensions of inverse length. The fluidised particles concentration is a dynamical 
variable in Vardoulakis et al. [3]. In subsequent works [12, 13, 21], the analysis is simplified 
treating this quantity as a constant, which is absorbed in the proportionality constant of the law and 
eliminated from the formulas. More recent works modify this law [30-33] setting the erosion rate 
proportional to the pressure gradient minus a cut-off value, in order to incorporate the 
phenomenology of the onset of sand production more directly into the modelling. In this work, the 
choice of the erosion model serves mostly illustration purposes. As discussed in the introduction, we 
shall follow Papamichos et al. [13] fairly closely, adopting the basic hydro-mechanical 
characteristics of the erosion model used there. 
The erosion model equation we shall use reads 


m 


=Ady(r,t)(1— pfa (7) 
Psolid 


where m is the rate of eroded mass production per unit volume, psojia is the mass density of the solid 
matrix of the porous material and A is a quantity with dimensions of inverse length and represents the 
strength of the erosive processes that lead to sand production. Parameter À is usually called sand 
production coefficient, but we shall prefer the name erosion strength coefficient, or merely erosion 
strength. Parameter J is an exponent-parameter of the model. Parameter 4,(r,f) is the erosion profile 
function defined as follows. 

The material is also subjected to fixed radial symmetric stress conditions at its inner and outer radius. 
Because of the imposed stress, plastic deformations are created around the hole. Let A be the depth of 
the plastic region. The /,(7,1) can be any smoothed step function such that equal to 1 (unity) within the 
plastic zone and 0 (zero) outside. Explicitly, we require 


Tout 
Jp (r, dr = ACi) (8) 


We choose an exponential function with these properties 


Ap(r,t) = exp[—A(t) “ [E0 + 1/a)]*x (r = rin)“] (9) 


(I(x) is the usual Euler gamma function.) The profile function /,(7,t) approaches a step function as a 
increases. The exponent a is fixed at the beginning of the analysis, and it is not tuned further as part 
of the calibration procedure. In this work, we shall use a=2 throughout. The shape of the profile 
function for various values of a is shown in Figure 2. 

The solid material mass conservation reads m = /p,,,g0@/Ot. Therefore, Eqn (7) takes the 
explicit form 


Copyright © 2015 John Wiley & Sons, Ltd. Int. J. Numer. Anal. Meth. Geomech. 2015; 39:2017-—2036 
DOI: 10.1002/nag 


aSU SUOWUOD ALAJ afquordde ayy Aq pausaaos are sapone YO ‘asn Jo sans 104 Areagry] IUUQ ÁM UO (suonrpuoo-pue-swa/wW03' AaytA\-AresqIauljuo//:sdyy) suoMIpUOD pue SUJAL AP 2ag “[EZOZ/OI/TE] UO Aresqry] ouuo Aart ‘Aresqr] uono JO Ausan Aq EYET SLU/ZOOI OI /OPfuoo Aaya Aresquaurpuo//:sdyy Woy papeojuMod ‘81 ‘STOT ‘ES869601 


2022 E. GRAVANIS, E. SARRIS AND P. PAPANASTASIOU 


1.0 


profile function 4, 
oO 
uw 


S 
(e) 


1 
(r—r,)/A 


Figure 2. The erosion profile function /, for a=2, 4, 8, 20, 100. 


e = 2% (r,t) (1 — fq (10) 


As implied by Eqn (8), the plastic region depth A changes with time as the stresses change with 
drawdown and depletion time. The rule under which A evolves will be described in the succeeding 
texts. The model is completed by introducing a prescribed dependence of the erosion strength 2 on 
A, also given later on. One should therefore bear in mind that 2 depends on time. 


2.3. Solution of the radial erosion flow 


The problem can be treated either under fixed pressure drop Ap or fixed flow rate Q. For the cases, we 
shall consider that the two conditions are practically equivalent. We shall work with the mathematically 
more interesting case of fixed pressure drop. In Section 6, we shall discuss the constant rate 
condition in the context of generalising the erosion model (Eqn (7)). We define the mathematical 
time T by 


dT = -A(nd 11 
= q Ad (11) 


where Ar=Prout— Fin. Then the erosion model (Eqn (7)) becomes 


õp Ar 


f 
A T)(1— @) (12) 


Its solution reads 


1 = 
#(r,T) = Po = gc (13a) 
(1+ 6- D0 - do" loat nar) 
Ar (T 
o(r,T) = 1 — (1 — Ø)exp -Z fo%lr, rar] (13b) 


for B#1 and f=1 respectively. Now, everything can be calculated in terms of d= ¢(r,T). A(T) is 
calculated from this function by Eqn (6). The corresponding physical time is calculated by 
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T aT 
t(T) = Arjo —> 14 
(r) = arbo aH (14) 
Then the eroded mass rate is given by 
: P Tou 
M = [mdV = pyig2eHA(Al, dp(r,t)(1 — p)ar (15) 


2.4. Calculation of the plastic region depth: erosion poro-elastoplastic solutions 


Once #(r,T) is known, one may integrate Darcy’s law (Eqn (2)) to obtain the evolution of the pressure 
profile. One finds 


r (1—¢)'d 
P= Pin + mh £ a = Pa EAT), ea (16) 


The first expression in Eqn (16), that is the one involving pin, can be effectively replaced by a 
polynomial approximation within the erosion region; this merely requires the numerical evaluation 
of the integral at two or three points. Outside the erosion region, the porosity remains essentially 
fixed at its initial value; hence, the second expression in Eqn (16) can be evaluated explicitly. (Eqn (16) 
should be understood according to these simplifying observations.) Thus, all calculations involving 
pressure can be performed explicitly in what follows. Throughout the region Fin<r<řrfouw, the stress 
components o,, Gg satisfy the equilibrium equation: 


d r r= 
r To (17) 
dr r 
What follows should be understood in the quasi-static sense. Therefore, explicit time dependence 
will be left understood. The whole region is divided into a plastic region, rin<r<rin+A, and an 
elastic region, ri,» + A<r</rou (Figure 1). In the elastic region, we have 


Oy = Op + (ALame F 2G)e, F ÅLame£0 (18a) 


09 = Op T ALame€r + (ALame F 2G)eo (18b) 
where the Lame parameters are given by 


E vG 
ae a ALame = 7—>~ 
5 S (1 — 2) 


2(1 +) (13) 


E is the Young modulus [MPa] and v is the Poisson ratio [—]. The strain components (radial strain €, 
and tangential strain ¢g) are expressed in terms of the radial displacement u by 


_ du 
Or 


Er 


yas (20) 
r 
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Then, for any given pressure profile p(r), the equilibrium (Eqn (17)) gives 


C2 a A ' 1 , 
u = Cyr + — — ———_ l], rplr )dr 21 
EG jae P(r) 
Cı and C, are integration constants, which are to be determined in terms of the input data. The 
parameter a is the Biot coefficient, and in our explicit calculations here, we shall set a= 1, meaning that 
we will neglect any compressibility effects. Employing the boundary condition G, etastic(“out) = Cout 
eliminates the constant C4, and one finds the expressions for the radial and tangential stresses: 


1 1 a2G L wipes 
r = 9ou C2 2G d 22 
ý Paepe = a T Alame + 2G > bad P(r ) í | i j 
1 1 a2G L Ga ee aa 
= Gon + C2 2G d 22b 
79 Pout + 2 | am F a a Atame + 2G po r2 frar pl ) r! i i 


For the creation of the plastic region, we have considered the Mohr—Coulomb failure criterion: 


cg — ap = So + K(a; — ap) (23a) 
a/n @® 
K = tan’ |7 +5 So = 2CVK (23b) 


where C is material cohesion [MPa], ® is the material friction angle [°] and So is the uniaxial 
compressive strength [MPa]. Equilibrium (Eqn (17)) and the boundary condition o(rin) = in give 


rK! S r£! ere ast 
Or = Fin RI Z|! 5 a(K nf“ fa (r) “p(r)dr (24) 


At this stage, there are two continuity conditions and two unknowns, the integration constant C2 and 
the location of the plastic zone boundary location R. The condition G,etastic(R) = Grplastic(R) gives C2 
explicitly in terms of R 


1 1 RK! So RX! 
czo- | = ont I! — Sein 


out R in Tin (25) 
oT ee ee ee a2G DE ae ee 
—a(K —1) r e (r ) p(r jar ~ hn FG erat PC jar! 
and 
O9 elastic (R) = OQplastic (R) (26) 


is solved numerically for R. Given the continuity of o, and p, the continuity of og follows from 
applying the Mohr-Coulomb failure criterion (Eqn (23a)) at r=R also on the elastic solutions 
(Eqns (22a) and (22b)). Note that Eqn (26) is purely algebraic. Thus, the plastic region depth is 
determined, A =R — rin. 


2.5. Degradation: Running A and the quasi-steady state 
The plastic region depth A is affected indirectly by the hydrodynamic evolution of the pressure profile, 


given in Eqn (16), which in turn affects the calculation of A, discussed in the preceding texts. This 
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effect is actually small. We now come to the law under which A becomes dynamical in a direct manner 
because of the deterioration of the material caused by erosion. 
By definition of the erosion profile function, we have that at any instant of time: 


Tout 
Ary ra 4r, T)dr = 1 (27) 
which is nothing but re-writing (Eqn (8)). We use this fact to define a spatial average of the porosity: 


BT) = A(T C,T) (r, T)dr 28) 


We introduce the running Young modulus and cohesion coefficients defined by 


1- (7) _ p1- 47) 
1— Øo ANRE 


expressing the degradation of the material. We assume that the friction is left unaffected by 
degradation. At the external confinement used in our work and Papamichos et al. [13], friction is 
possible to be mobilised and increase as a function of the plastic shear strain. However, plastic shear 
strains also cause deformation on the grains, causing the creation of cracks that directly reflects in a 
cohesion-softening behaviour. As the cementation of the natural or the synthetic sandstones 
gradually decreases, the cohesive strength of the material weakens. In such cases, micro-cracks 
develop first (which is usually the case) and it is reasonable to assume a cohesion-softening 
behaviour without the friction hardening. The later assumption was also considered in Stavropoulou 
et al. [21]. 

A(T) is calculated through Eqns (16-26) using the running values of the material parameters. 
Explicitly, the system is imagined to evolve through a succession of quasi-static time-steps. In the 
duration of each step, the value of A as well as of E and C are constant and serve to determine the 


E(T) =E (29) 


value of A associated with the next time step. Simultaneously (7,7) and therefore (T) are 
calculated in each time step by Eqn (13), taking into the values of A of the all previous time steps, 
as we explain in Section 2.7. This determines also the values of E and C associated with the next 
time step. This rule makes the plastic region depth dynamical in a direct manner and also couples 
explicitly, and strongly, the hydrodynamic part of the erosion model to its material component. 

One should note that the values for the Young modulus and cohesion calculated by Eqn (29) are 
rather nominal values through which one models the degradation of the material in order to estimate 
sand production. Subsequently, the derived stress and pressure fields as well as the plastic region 
depth are similarly approximate. These variables are rather part of an all too approximate description 
of erosion, starting with the model in Eqn (6) or any similar model, operating as auxiliary quantities, 
which provide a very rough description of the actual local physics of the erosion process. 


2.6. Switching on the erosion strength i 


The erosion strength is the less understood part of models that involve equations similar to Eqn (7); 
this applies to a large proportion of the existing sand production models. (For sand production models 
based on different principles, more intimately related to the mechanics of the material, refer to the 
discussion in the review work of Rahmati et al. [29].) The reason why / is poorly understood is that 
the associated models are primarily hydro-dynamical in nature. In fact, they are purely hydro- 
dynamical when 2 is an absolute constant. The mechanics of the material is usually introduced by 
turning the erosion strength / to a function A(é,) of the local equivalent plastic strain €, (e.g. refer to 
Rahmati ef al. [29]). This method, introduced in Papamichos and Stavropoulou and Papamichos 
et al. [12, 13] induces association of the erosion region with the plastic region, thereby localising 
erosion. Simultaneously, it causes a gradual switching on of the erosion strength, modulated by the 
plasticity. 

In the present work, localisation is effected by the profile function 4,(r,1), which is controlled by the 
depth of the plastic region. The magnitude of the erosion strength J is still a constant. As we shall 
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discuss in Section 5, A is increasing during the initial stages of the erosion process, then reaches a stage 
where it can be approximated as constant, and eventually reaches a stage where it is decreasing. 
Therefore, the time variation of the 2 must be introduced independently. Presumably, this latter stage 
is the most difficult to understand as it seems to be related to modes of failure of the material lying 
beyond the modelling of the material presented here. It is not within the scope of this work to 
improve the mechanics concerning the material aspect of the sand production models. Thus, and for 
illustration purposes only, we shall model the initial to middle stage variation of 2 by introducing 
time dependence to 4, requiring that it depends on the running plastic region depth A by the 
following simple rule: 


A-A 
2o + (A2 — 20) i ASA, 


A(A) = Ay — Ao’ (30) 
Ad, A 2 Ay 


The numbers 20, 42 and A, are model parameters that must be tuned by calibration. Figure 3 shows 
the graphical representation of this relation. 

In words, Eqn (30) expresses the following. The erosion strength starts at an initial value 2p and 
grows linearly with A. If A exceeds a (threshold) value A,, then the strength 2 becomes constant 
and equal to 22 (reaches a plateau). The threshold value A, is quantified in terms of the initial depth 
Ao, which is calculated from the equilibrium equations using initial values of the input data (refer to 
also the next section). Therefore, À is essentially a function of the ratio A/Ag. 


2.7. Solution procedure 


The formulas given in the preceding texts reduce the derivation of explicit results to the calculation of 
certain integrals and the solution of an algebraic equation through a series of time steps. The 
independent time variable is the mathematical time T. T is discretised according to a suitable time 
step OT; that is, we set T=idT, where i is an integer. 

At T=0, the porosity field is everywhere known, @= ġo. Therefore, the quantity A from Eqn (6) and 
the pressure field from Eqn (16) are calculated easily. Then the equilibrium equations of Section 2.4 
allow the determination of the initial plastic region depth Ag= A(T=0) by solving Eqn (26). 

At time T=6T, the porosity field is calculated from Eqn (13a). The formula involves the integral 


{pant T)dT, which is approximated by 4,(r,0) ôT. The profile factor 2,(r,T) is given by Eqn (9). 
The 4,(r,0) is then known, as Ag=A(T=0) is known from the previous time step. The quantity A 
and the pressure field are calculated from Eqns (6) and (16), and the equilibrium equations provide 
the new value of the plastic region depth A; = A(T=6T) by Eqn (26). For the latter calculation, the 
initial values of the Young modulus and the cohesion are used. Equations (28) and (29) allow then 
the calculation of the running Young modulus and cohesion, E; = E(T=0T) and C;=C(T=6T), to be 
used in the next time step along with A. 


À 


Figure 3. Graphical representation of Eqn (30). 
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At the general time T=iôT, the porosity field is calculated from Eqn (13a) once the integral 


{oap(r, T)dT is approximated by {4,(r,0)+...+ 4,(7,@— 1)6T)} OT. This involves the values of the depth 
A at all the previous time steps, which are known. As explained in the preceding texts, the quantity A, 
the pressure field and the new value of the plastic region depth A;= A(7=i67) can then be calculated. 
By Eqns (28) and (29), the Young modulus and cohesion, E;=E(T=idT) and C;=C(T=i0T), to be used 
in the next time step along with A,. 

Given the information collected as explained, Eqns (14) and (15) allow the translation to real time t 
and the calculation of sand production for all time steps. Also, the radial and tangential stresses g, and 
og can be calculated from the formulas of Section 2.4. Presumably, the discretised version of the model 
discussed in the preceding texts could be very well regarded as the actual proposed model in this work. 


3. COMPARISONS BETWEEN MODEL PREDICTION AND EXPERIMENTS 


In this section, we present the results from the analysis conducted for the modelling of the sand 
production in the hollow cylinder test problem. The rock is treated, in accordance with our 
formalism, as an elastic-plastic Mohr—Coulomb solid infiltrated with fluid which obeys the Darcy law. 

The mechanical properties of the synthetic sandstone used in the hollow cylinder test of Papamichos 
et al. [13] were studied experimentally in Tronvoll et al. [2] for external stresses up to 3.5 MPa. The 
test presented in Papamichos ef al. [13] involved external stresses above 7.5MPa; for the 
simulations, the authors used a non-linear elastic-plastic model [30] calibrated for synthetic 
sandstone according to the data of Tronvoll et al. [2]. Here we shall use the mechanical properties 
of the Red Wildmoor sandstone calibrated for the Mohr—Coulomb yield criterion after the work of 
Sulem et al. [34]. Any difference between these material data and the actual data of the synthetic 
sandstone at the stress levels of interest will not affect qualitatively the issues of primary concern in 
this work; that is, whether the basic physics underlying the models in Papamichos and Stavropoulou 
and Papamichos et al. [12, 13] is maintained in the present model. Nonetheless, in Section 6, we 
shall return to investigate parametrically the effect of the cohesion C and friction angle ®, which 
predominantly affect plasticity and therefore the sand production, through back analysis. The input 
and model parameters upon which the computations were performed are given in Table I. 


Table I. Model properties and input data. 


Variable Value 


Geometric properties 


Hollow cylinder internal radius, 7;, [m] 0.01 

Hollow cylinder external radius, Fout [m] 0.1 

Cylinder height, H [m] 0.2 
Porous rock and fluid properties 

Young modulus, £ [MPa] 6750 

Poisson ratio, v [—] 0.19 

Cohesion, C [MPa] 3.7 

Friction angle, ® [°] 37.4 

Initial porosity, do [—] 0.3 

Initial rock permeability, k [md] 500 

Solids density, psoria [kg/m?] 2640 

Kozeny—Carman parameter, k, [m7] 8.956E-12 

Dynamic viscosity, u [MPas] 5.0E-9 

Model Parameters 

Erosion strength, 2 [m7] Eqn (30) 

Initial erosion strength, 2o [mt] 0.42 Az 

Maximum erosion strength, 22 [m] 0.088 

Initial plastic region depth, ^o [m] equilibrium equations 

Threshold depth of the À plateau, A, [m] 1.1 Ao 


Model exponent in Eqn (7), 2 [—] 


1 
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Based on the test conditions presented by Papamichos et al. [13], we have performed simulations 
with different external radial stresses Cou =7.5, 8, 9, 10, 11 MPa, while the radial inner stress Gin 
was kept to zero, and to an external pressure of Pow =0.15 MPa (with internal pore pressure pin 
equal to zero), which corresponds to a flow rate Q=0.5 l/min. 

In Figure 4, the sand production curves produced by the model under the fixed pressure drop 
Ap=0.15 MPa for all the different values of external radial stress are compared with the experimental 
data of Papamichos et al. [13]. All the physical parameters are taken from Papamichos et al. [13], 
while the model parameters, associated with the erosion strength 2 (shown in Table I) are tuned so 
that to obtain a reasonable fit with the experimental data. 

Most of the tuning is required in order to capture the initial stages of the process. The tuning might 
be slightly excessive as the data may be distorted by the gradual change of flow rate in the experiment 
of Papamichos et al. [13]. Nonetheless, we shall take the data presented in Papamichos et al. [13] at 
face value and attempt to calibrate according to them. No dependence of 4 on the external stress is 
required: Once, the model was calibrated for a certain reference value of external stress (11 MPa), 
that is the predicted curve to fit reasonably well with the experimental points, thet rest of predicted 
curves fall in place. One should note that at 7.5 MPa, sand production declines strongly, as if there 
were an external stress cut off value for the initiation of sand production. One should observe this 
aspect of the experimental data shown in Figure 4 in conjunction with the progression of the plastic 
region depth A, as predicted by the model for the same cases, shown in Figure 5. One observes that 
especially at the very initial stage, the plastic region depth for the case 7.5 MPa differs by more than 
an order of magnitude from the depth of the next case (8 MPa), while the other four cases lie within 
an order of magnitude. The initial values A are mostly a result of the standard equilibrium equations 
combined with the Mohr—Coulomb criterion than of any of them are components of the model. We 
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Figure 4. Comparison between the curves predicted by the model (denoted as ‘M’ - continuous lines) with the 
test data of Papamichos et al. [13] (denoted as ‘L’ — marker points) for fixed pressure drop Ap=0.15 MPa 
(flow rate Q =0.5 l/min). 
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Figure 5. Progression of the plastic region depth A as predicted by the model for the cases of 7.5, 8, 9, 10 
and 11-MPa external radial stresses. 
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consider this observation as supporting the empirical basis of the localisation method of Papamichos 
and Stavropoulou and Papamichos et al. [12, 13], which we have also adopted here: Plasticity 
controls erosion. In regard with the late time stages of the process, one should note that the 
experimental curves in Figure 4 exhibit break points, that is parts of the curve where slope has 
decreased suddenly. This phase of the process is associated with material failure not covered by the 
Mohr-—Coulomb criterion and therefore not described by the model. 

In all, these results are in accordance with the findings of Papamichos et al. [13], verifying that the 
present model carries essentially the same physics as the one used there. The present model may have a 
little more freedom, as the switching on of the erosion strength 4 with time, given by Eqn (30), is not 
modelled simultaneously with localisation, which is effected by the profile function in Eqn (7). The 
general statement is that hydrodynamic erosion models of the general type (Eqn (7)), which localise 
erosion through plasticity, perform reasonably well under different stress levels, although they are 
compromised because of possible additional modes of material failure during the process. 

In Figure 6, we present the complete solution of the coupled hydro-mechanical erosion problem. 
The solution includes (a) the prediction of the cumulative sand production, the hydro-mechanical 
response with time of the physical fields of (b) porosity, (c) pore pressure and (d) fluid (Darcy) flux 
q; the mechanical response of the hollow cylinder (e) radial stress and (f) tangential stress, for the 
specific case of external radial stress at 11 MPa with pressure drop Ap=0.15MPa. The curves 
correspond to different stages of the erosion process differing by 1000s. The final erosion stage is 
set at 9000s. 
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Figure 6. Sand production hydro-mechanical solution according to the proposed model for external stress 

level 11 MPa and pressure drop 0.15 MPa: (a) cumulative sand production, (b) porosity profiles, (c) pore 

pressure profiles, (d) fluid flux profiles along the hollow cylinder, (e) radial stress profiles and (f) tangential 
stress profiles from time 0 to 9000s for every 1000-s interval. 
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The effect of erosion is primarily encoded in the variation of the porosity field. Porosity starts from 
the uniform initial value 0.3 and increases near the inner hole. This causes material softening, which in 
turn causes the plastic region depth A to progress and the porosity increase advances deeper into the 
rock. The flattening of pressure curves observed near the hole, that is the decrease of the local 
pressure drop in the vicinity of the inner radius, is of course a direct effect of erosion as encoded in 
the porosity field (permeability increases with porosity by the Kozeny—Carman equation). Softening 
of the material, along with the advancement of the plastic zone, is reflected in the tangential stress 
solution, where stresses near the inner radius decrease significantly as a result of the coupling 
between hydrodynamics, plastic yielding and material degradation (softening). The effect of these 
three mechanisms coupled together can be described as follows. Hydrodynamics and degradation 
erode and soften the material, and plastic yielding provides that more material is made available for 
erosion and softening. 

The behaviour of the erosion model (Eqn (7)) when varying the pressure drop Ap, or equivalently 
the flow rate Q, is summarised as follows. In Eqn (7), the local mass production is proportional to 
the Darcy flux. This induces a linear dependence between the total sand production and the flow rate 
Q. Indeed, we may introduce a time scale out of the Darcy flux (velocity) scale (k,/u)(Ap/Ar) and 
the erosion strength that has dimensions of inverse length: 


“Ar 
= 1 
T= IRAD a 


The parameter 4 in this formula may be regarded as the maximum value 42. Figure 7 shows the 
sand production curves derived by the model for different values of pressure drop, Ap=0.15.,..., 
0.90 MPa with a constant external stress level of 11 MPa against the scaled time t/t. From Figure 7, 
one observes that all production curves are essentially identical. The difference is mainly because of 
the evolving A during the initial stage of the process. Therefore, the sand production is indeed 
proportional to the flow rate under the Eqn (7); that is, if we were to calibrate the model for a certain 
reference value of Q, then we can predict sand production for other values of the flow rate. In 
Papamichos et al. [13—28], it is argued that this proportionality is supported (very) approximately by 
experiment. Nonetheless, considering the data of Papamichos et al. [13] in scaled time t/t suggests 
that this is too approximate and possibly not entirely correct, in the sense that different Q curves 
(for the same external stress) do not visibly collapse to each other at an acceptable level. At any rate, 
we shall regard the g dependence of the erosion models as an open issue. We will revisit this matter 
in Section 6. 
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Figure 7. Predicted sand production curves for external stress 11 MPa and the different pressure drop level 
(different flow rate) in scaled time. The curves are practically self-similar. 
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4. INFLUENCE OF MODEL PARAMETERS ON SAND PREDICTION 


The model parameters are internal quantities to the model without an established physical significance. 
As mentioned in Section 2.6, these quantities are the initial erosion strength Ao, the maximum erosion 
strength 22 and the threshold depth A, of Eqn (30). Furthermore, it is useful to add to them the 
exponent p of the model (Eqn (7)) that changes the porosity of the material. 

It is instructive to study the effect of these parameters on the sand production curves, that is how the 
shape of the curve is changed by varying the model parameters. Thus, we perturb the model parameters 
around the values obtained after calibrating with the case of external radial stress 11 MPa and pressure 
drop 0.15 MPa with the test data of Papamichos et al. [13]. The results from the analysis of the model 
parameters are shown in Figure 8. For each investigation, the model sand production curves are 
compared with the corresponding experimental curve. 

The effect of the maximum erosion strength 42 (Figure 8a) and the effect of exponent / (Figure 8d) 
are to change the overall pace of sand production in the post-initial stage of the process, while the 
exponent f affects also slightly the curving. Figure 8b shows that the threshold ^, mostly delays the 
post-initial stages, and the predicted sand production curves behave in the post-initial stages as 
curves of smaller / than 42. Finally, the initial strength 29, or better the fraction Ap/A2, has the same 
effect with the threshold A,, parameter along with a suppression of the initial sand production. 


5. EROSION STRENGTH 2 FROM EXPERIMENTAL DATA 


All our investigations are based on the rule governing the evolution of the erosion strength À during the 
erosion process, given by Eqn (30). This is an empirical and crude way to simulate the behaviour of 4 
during erosion, so the model sand production curves to imitate reasonably the experimentally deduced 
production curves under given condition (as it is shown in Figure 4). Then, the following interesting 
question arises naturally. How å need to vary over time, that is how its exact pattern looks like 
during the erosion process in order to produce the experimental sand production curves as model 
production curves? 

The answer to this question can be given by means of back analysis of the experimental data. The 
experimental curve (i.e. the set of points) of each case studied in Section 3 may be regarded as a 
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Figure 8. Influence of model parameters on the shape of the sand production curve with external stress level 
11 MPa and pressure drop 0.15 MPa: (a) maximum erosion strength 22, (b) threshold depth A,, (c) ratio of 
initial to maximum erosion strength 2/42 and (d) model exponent coefficient /. 
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piecewise constant sand production rate M curve, by connecting the experimental points with straight 
lines. Then, at each time step, A can be deduced from Eqn (15), which is now used in the inverse way— 
from a known sand production curve M, we deduce A over time. With a value of 4 at hand for each time 
instant, one may proceed as discussed in Section 2.7. 

In Figure 9, we plot the erosion strength 4 and the dimensionless product AA (erosion 
strength x plastic region depth) over time, for the five cases considered in Section 3. If erosion can 
indeed be associated with plasticity, as it is the physical basis of erosion localisation in the type of 
models we use here, then 2A, which is a pure number, must necessarily be an important quantity. In 
the present analysis, we attempt to reveal aspects of its behaviour, in the hope that they will provide 
useful information for shaping an understanding of this quantity in the future. Presumably, when the 
exponent f of the model in Eqn (7) is set to 1, the sand production rate formula (Eqn (15)) can be 
written in the form 


M = (AA) Psia (1 — $) (32) 


using the spatial average ¢ of the porosity field within the erosion region, defined in Eqn (28). Thus, the 
number AA operates as the natural erosion strength for the total sand production rate, while Eqn (7) 
involves the local production rate. 

Inspecting the curves of the Figures 9a and 9b, one observes that both quantities, the A and AA 
graphs exhibit an initial stage where the erosion strength increases and a final stage where the 
erosion strength decreases, at least as a general tendency. As mentioned in Section 2.6, the erosion 
strength Å is forced to partially carry the burden of the mechanics of the material. One should note 
that although there are strong fluctuations, the values of A during the intermediate stage congregate 
somewhere below 0.1m~!. This is consistent with the intermediate stage value 1,=0.088 m~! 
(Table I), which we have used in order to obtain a reasonably good fit of the experimental curves. 
The AA curves cannot converge in the intermediate stage, because the depth A is different for each 
stress level case. Note that the excessive initial fluctuations of A, especially for the lower values of 
the external stress, are explained by the small initial plastic depth A in those cases. Therefore, by 
forming the product AA, those fluctuations are smoothed out, providing a more faithful description 
of the erosion strength in the initial stages than / itself. Presumably, any model, which localises 
erosion through plasticity, is bound to behave in a similar way. 

The previous comments are qualitative in nature. The most immediate aim of the back analysis is to 
reveal a general behaviour for A, specific enough from the quantitative point of view, so that to be able 
to improve the postulated rule 2(A) of Section 2.6, or simply replace it with something more 
appropriate. The difficulty with this idea is that fluctuations are large and the intervals of A for each 
stress level differ significantly and there is no reliable way to obtain a working approximation for 
the evolution of 4 from these data. 

On the other hand, back analysis may be used to investigate the effect of varying the input 
parameters on the estimation of 2 when the values of these parameters are poorly known or poorly 
estimated. As discussed in Section 3, this is the case of the parameters we have used in this work, 
which describe the mechanical properties of the synthetic sandstone of the hollow cylinder test of 
Papamichos et al. [13]. Thus, we calculate the mean value of 4 (over its course of evolution) for the 
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Figure 9. Behaviour of (a) erosion strength and (b) product of erosion strength with plastic zone depth with 
respect to time. 
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cases Cou =7.5 and 11 MPa (minimum and maximum external stress) for a range of values of the 
cohesion C and friction angle ®. The result is shown in Figures 10a and 10b for the minimum and 
maximum value of external stress respectively. We have considered values from 15 to 50° for the 
friction angle and values from 3 to 5MPa for the cohesion. The curves shown in Figure 10 are 
labelled by cohesion values that differ by 0.2 MPa. 

The larger the values of the cohesion C and the friction angle ® are the more the material resists to 
plasticity. Hence, as these values increase, the plastic region depth decreases, leading to a decrease in 
sand production. Consulting Figure 5, we may recall that the initial value of the plastic region depth Ap 
for the case Gout = 7.5 MPa is much lower than the value of Ag for the cases of Cout = 8 MPa and higher 
(for the C and ® values quoted in Table I). This is consistent with the findings of Papamichos et al. [13] 
that at 7.5-MPa external stress, the sand production drops significantly (as shown by the data in Figure 1). 
It is therefore reasonable to continue the curves in Figure 10 until the initial plastic region depth drops 
below a certain reference value, which is regarded as the minimum plasticity for sand production 
(which we set rather arbitrarily to 5x 10~°m, consulting again Figure 5). Based on these remarks, the 
first result obtained from Figure 10a is that the values of C and ® along the notional line formed by the 
endpoints of the curves correspond to very small sand production. Hence, we conclude that the actual 
values of the cohesion and friction angle must be along, or in the vicinity of that notional line. 
Secondly, one observes that for the entire range of both parameters along this line, the mean value of 4 
lies roughly between 0.07 and 0.105m7'; that is, there is no significant deviation from the value 
0.088 m~' (for the maximum erosion strength 22) we found in the preceding texts through calibration 
for the C and ® values quoted in Table I. 

In Figure 10b, there is a sequence of points where the curves of same C for the two different external 
stress levels cross each other. At such points, the same value of mean 4 is obtained for both external 
stress levels when calibrating with the same values of C and ®, as we have performed previously 
with the values quoted in Table I. It is almost needless to say that the same, or at least, roughly the 
same value of mean 4 should arise for different stress levels if an erosion model is to be useful at 
all. The mean / values obtained this way range from 0.025 to 0.09m!. This is a relatively large 
range, suggesting that one must calibrate with a different set of C and ® values for the different 
external stress levels. However, by the same token, if one insists on a A value between 0.07 and 
0.105m~', which seems to be suggested by the analysis of cou =7.5 MPa case, there is a fairly large 
range of C and © values consistent with those values of 2 in the oo4.=11 MPa case. The actual 
values of C and ® for the 11-MPa external stress may very well be within that range, most possibly 
towards the large C end. 


6. GENERALISATIONS OF THE EROSION MODEL AND EXTENSION OF THE WORK 


We have studied the erosion problem using the model given by Eqn (7) under fixed pore pressure 
difference between the inner and outer radius of the hollow cylinder test. As mentioned in Section 
2.3, the problem could have been equally well solved under fixed radial flow rate Q. Adopting the 
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Figure 10. Mean value of erosion strength 4 for the cases (a) cout = 7.5 MPa and (b) 11 MPa. The curves are 
labelled by the values of the cohesion C in range 3-5 MPa differing by 0.2 MPa. The curves of the case (a) 
are also plotted in dashed line for comparison with the case (b). 
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fixed Q condition, the solution given by Eqn (13) can be written directly in the physical time ¢. This, in 
turn, allows treating more general erosion models. 
Let us generalise the model (Eqn (7)) as follows: 


"= diy (r, (1 AFF) (33) 

Psolid 
where f(q) is a suitably chosen function of the Darcy flux q. For example, we may set f(q)= 4 — qer for 
q> er Otherwise 0, with qer being a cut-off for the flux densities under which erosion may occur. 
Clearly, qer is an additional model parameter. Models of this kind have been suggested and 
employed in the simulations of sand production in the recent works of Detournay et al., Detournay, 
Papamichos, and Azadbakht et al. [30-33]. For constant flow rate Q, the Darcy flux (refer to Eqn (5) 
in Section 2.3) 


Q 
= 34 
7 2aHr 83 
is independent of time. Then the counterpart of solution (Eqn (13a)) in the present model reads 
ie 
(r,t) =1 di A 176-1) 39) 
(1+ 6- D0 ~ y fr Aarna) 
The production rate formula (Eqn (15)) also becomes 
; ‘ie B Q 
M = Pooig2tH J, Ap(r,t)(1 — AVF ET rdr (36) 


Only obvious modifications are required in the formalism of Section 2 in order conform to the new 
relations (Eqns (33-—36)). The interesting attributes of the model given by Eqn (33) can be described as 
follows. 

In Section 3, we discussed the linearity of the model given by Eqn (7) with respect to the Darcy flux q 
as a possibly weak characteristic of the model. Detailed work is then required to show that the 
generalised model (Eqn (33)), for any specific choice of function f, performs well. That is, one should 
examine whether the predicted sand production agrees well with the experiments for different flow 
rates Q, while maintaining the erosion strength 4 and the additional model parameters, such as qer, 
independent of Q. At any rate, the novelty of the model (Eqn (33)), that is the non-trivial dependence 
on q, seems to be the right way to proceed. We shall defer the analysis of such models to future work. 


7. SUMMARY AND COMMENTS 


In the present study, we analyse the problem of sand production. We introduce an erosion model based 
on coupling hydrodynamic erosion with mechanical deformation of the rock under external stress. The 
two processes are linked by softening of the material; that is, the elastic stiffness and the cohesion of 
the rock are modelled to decrease with increasing porosity. In this description, sand production 
initiates as the material near the surface of the inner cavity fails mechanically by the applied stress 
field. The failed material erodes away while softening progresses deeper into the rock making more 
material available for erosion. The introduced hydro-mechanical erosion model while maintaining 
the physics usually incorporated in erosion models presents the advantage that it is more susceptible 
to detailed mathematical analysis. As a result, a mathematical framework is set up from which the 
explicit results can be deduced with relative ease, that is without the need to employ advanced 
numerical methods or use specialised software. 
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For calibration and validation purposes, we set the theoretically derived sand production curves 
against experimental data. The performance of the model rests on the dependence of the erosion 
strength J on the plastic zone depth A, discussed in Section 2.6. The dependence of A on plasticity 
parameters is not necessary for the hydro-mechanical coupling in the present model; it only serves 
of additional flexibility. The function 1(A) involves three free parameters that are to be tuned by 
calibration. The results show reasonable agreement between predictions and experiment. This 
analysis verifies the claim that the basic physical ideas usually incorporated in hydro-mechanical 
models have been successfully embodied in the present model. 

The dependence of the erosion strength 2 on plasticity parameters, in our case the plastic zone A, 
carries information about the effect of the mechanical deformations on erosion not covered by the 
hydro-mechanical coupling. This is why A must be a dynamical variable and not an absolute 
constant. The variation of 2 with time can be read off from the experimental data by back analysis. 
Performing this analysis, we find that A exhibits a tendency to increases during the early stages of 
erosion and decrease during the late stages. Although the experimental noise is large and no definite 
general behaviour of À can be deduced from the results of this analysis, the results themselves can 
be practically useful as tabulated information about the material under a given set of conditions. 

One particular weakness about the type of the hydro-mechanical sand production models we 
consider here is their behaviour under different fluid flow rate conditions, that is that they behave 
linearly: The sand production is practically proportional to the fluid flow rate, which—in the best 
case—is a crude approximation of actual sand production. We suggest generalisations of the model 
that may have sufficient structure to deal with this difficulty and point out that these more general 
models can be equally well accommodated in the constructed mathematical framework. Detailed 
examination of their performance is left for future work. 
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